function [sgb, ...
          equalorder, ...
          Ugridreduced, sgreduced,...
          id_Ugridreduced,...
          s_id_gridreduced, sU_id_gridreduced,...
          n_U_sets, Tcoord, indices_pairs,...
          numeqconstraints,numzeros,coordrowseq, d3, e5, f5, g3, h5, i5, l3,m5, n5,  fillAeq]=useful_anyparam3(u0grid, u1grid, u2grid, sg)     
%% U sets
n_U_sets=3; %number of U-sets (number of dimensions of the distribution of the taste shock differences)
Ugrid=cell(1,n_U_sets); %cell 1xn_U_sets
Ugrid{1}=[u0grid-u1grid -u0grid+u1grid u0grid-u2grid -u0grid+u2grid u2grid-u1grid -u2grid+u1grid zeros(sg,1) Inf*ones(sg,1) -Inf*ones(sg,1)];
Ugrid{2}=[u0grid-u2grid -u0grid+u2grid u0grid-u1grid -u0grid+u1grid u2grid-u1grid -u2grid+u1grid zeros(sg,1) Inf*ones(sg,1) -Inf*ones(sg,1)]; 
Ugrid{3}=[Inf*ones(sg,1) -Inf*ones(sg,1) u2grid-u1grid -u2grid+u1grid u0grid-u1grid -u0grid+u1grid u0grid-u2grid -u0grid+u2grid zeros(sg,1)]; 

sgb=0;     

%% Reduce number of grid points to check by applying our partitioning arguments
%equalorder (sgx1) assigns the same index to parameters vectors such that Ugrid{1}, Ugrid{2}, ..., Ugrid{n_U_sets} have the same order

%%%%%%%% Criterion 1 %%%%%%%%
d=cell(n_U_sets,1);
index=cell(n_U_sets,1);
for j=1:n_U_sets 
    [sortedUjgrid,index{j}] = sort(Ugrid{j},2); 
    %sortedUjgrid gives Ugrid{j} with each row sorted from smallest to largest; it is sgxn_columns_Ujgrid
    %for each row i of sortedUjgrid and for k=1,2,...,n_columns_Ujgrid, index{j}(i,k) gives the column position in Ugrid{j}(i,:) of sortedUjgrid(i,k); it is sgxn_columns_Ujgrid
    d{j} = diff(sortedUjgrid,[],2) == 0; %sg x (n_columns_Ujgrid-1)
    %Suppose that Ugrid{j} has 3 columns. Then, for each row i of sortedUjgrid, d{j} will be: 
    %[1 1] if sortedUgridj(i,2)-sortedUgridj(i,1)=0 and sortedUgridj(i,3)-sortedUgridj(i,2)=0
    %[1 0] if sortedUgridj(i,2)-sortedUgridj(i,1)=0 and sortedUgridj(i,3)-sortedUgridj(i,2)~=0
    %[0 1] if sortedUgridj(i,2)-sortedUgridj(i,1)~=0 and sortedUgridj(i,3)-sortedUgridj(i,2)=0
    %[0 0] if sortedUgridj(i,2)-sortedUgridj(i,1)~=0 and sortedUgridj(i,3)-sortedUgridj(i,2)~=0
end
%%%%%%%% Criterion 2 %%%%%%%%
n_columns_Ujgrid=9;
Ugrid_sum=zeros(sg,n_columns_Ujgrid+n_columns_Ujgrid^2);
for g=1:sg
Ugrid_sum(g,:)=[Ugrid{1}(g,:) reshape((Ugrid{2}(g,:)+Ugrid{3}(g,:).')',[],1)'];
end

[sortedUgrid_sum,index_sum] = sort(Ugrid_sum,2); 
d_sum = diff(sortedUgrid_sum,[],2) == 0; %sg x (n_columns_Ujgrid*n_U_sets-1)

%%%%%%%% Combine %%%%%%%%
[~,~,equalorder] = unique([index_sum d_sum [index{:}] [d{:}] ],'rows', 'stable'); %sgx1 

%Reduced grid of parameter vectors and U sets
[~,index_final,~]=unique(equalorder, 'stable'); %one row from equalorder per ordering group together with its index
Ugridreduced=cell(1,n_U_sets);
for j=1:n_U_sets 
    Ugridreduced{j}=Ugrid{j}(index_final,:); %sgreduced x n_columns_Ujgrid
end
sgreduced=size(Ugridreduced{j},1); %size of the reduced grid of parameter values


%% Construct indices to keep track of equal elements inside Ugridreduced{j} for j=1,...,n_U_sets
id_Ugridreduced=cell(n_U_sets,1);
for j=1:n_U_sets 
    id_Ugridreduced{j} = NaN(size(Ugridreduced{j})); % %sgreduced x n_columns_Ujgrid
    for k = 1:sgreduced
        [~, ~, id_Ugridreduced{j}(k,:)] = unique(Ugridreduced{j}(k,:), 'stable'); 
    end
end


s_id_gridreduced=zeros(sgreduced, n_U_sets); %useful to set number of unknowns of the linear programming
for j=1:n_U_sets
    s_id_gridreduced(:,j)=max(id_Ugridreduced{j},[],2); %number of distinct elements of Ugridreduced{j}
end
sU_id_gridreduced=prod(s_id_gridreduced,2); %sgreducedx1 reporting the cardinality of the Cartesian product of the U sets
                                            %sU_id_gridreduced(g) is the number of unknowns of the LP for each g-th parameter value
                                            
s_id_gridreduced_2=ones(1, n_U_sets)*n_columns_Ujgrid; %number of elements of Ugridreduced{j} (count as separate even equal elements)
sU_id_gridreduced_2=prod(s_id_gridreduced_2); %cardinality of the Cartesian product of the U sets     
  

%% Invariant components of the linear programming problem 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% Inequality constraints %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

vectors = arrayfun(@(x) {1:x}, s_id_gridreduced_2); 
n = numel(vectors);
Tcoord_temp = cell(1,n);
[Tcoord_temp{:}] = ndgrid(vectors{:}); 
Tcoord_temp = cat(n+1, Tcoord_temp{:});
Tcoord = reshape(Tcoord_temp,[],n); %matrix of dimension sU_id_gridreduced_2xn_U_sets reporting the column coordinates of all the possible n_U_sets-tuples from the U-sets
%For example if n_U_sets=3 and s_id_gridreduced_2 is [2 3 1], then 
%[ca, cb, cc] = ndgrid(1:2, 1:3, 1:1);
%Tcoord = [ca(:), cb(:), cc(:)];   
%Tcoord= [1     1     1
%         2     1     1
%         1     2     1
%         2     2     1
%         1     3     1
%         2     3     1];
indices_pairs=pairIndices(sU_id_gridreduced_2); %row indices of all possible unordered pairs of rows from the matrix obtained by taking the Cartesian product of the U sets
                                         %[sU_id_gridreduced_2*(sU_id_gridreduced_2-1)/2]x[2]   
%Continuing the example above
%indices_pairs=[1 1; 1 2; ...; 1 27; 2 3; ...; 2 27; ... 26 27]


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% Equality constraints %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

numeqconstraints=3+(9+3*2*sgb+9+3*2*sgb-1)*(9+3*2*sgb-1)+(9+3*2*sgb)^2+1+(3+1+sgb*3)*3+(2*7)+(2*2*sgb)*3;  
numzeros=(9+3*2*sgb+9+3*2*sgb-1)*(9+3*2*sgb-1)+(9+3*2*sgb)^2;

coordrowseq=[1 1 1 ... %observational equivalence
             2 2   ...
             3     ...
             4:1:3+numzeros ... %zeros
             3+numzeros+1 ... %ones
             kron(3+numzeros+1+1:1:3+numzeros+1+(3+sgb*3)*3, ones(1,2)) ... %symmetry
             3+numzeros+1+(3+sgb*3)*3+1:1:3+numzeros+1+(3+sgb*3)*3+3 ...
             kron((3+numzeros+1+(3+sgb*3)*3+3+1:1:3+numzeros+1+(3+sgb*3)*3+3+(2*7)+(2*2*sgb)*3),ones(1,2))]; %identical  

         
%needed to construct coordcolumnseq    
%used to construct coordcolumnseq
d1=kron((10:1:9+2*sgb), ones(1,2));
d2=8*ones(1,2*sgb*2);
d3=[d1;d2];

e1=9+2*sgb+1:1:9+2*sgb+2*sgb;
e2=8*ones(1,2*sgb);
e3=[e1;e2];
e3=e3(:).';
e4=8*ones(1,2*sgb+2*sgb);
e5=[e4; e3];

f1=ones(1,2*sgb);
f2=9+2*sgb+1:1:9+2*sgb+2*sgb;
f3=[f1;f2];
f3=f3(:).';
f4=ones(1,2*sgb+2*sgb);
f5=[f4; f3];

g1=kron(9+2*sgb+1:1:9+2*sgb+2*sgb, ones(1,2));
g2=8*ones(1,2*sgb+2*sgb);
g3=[g1;g2];

h1=10:1:9+2*sgb;
h2=8*ones(1,2*sgb);
h3=[h1;h2];
h3=h3(:).';
h4=8*ones(1,2*sgb+2*sgb);
h5=[h4; h3];

i1=ones(1,2*sgb);
i2=9+4*sgb+1:1:9+4*sgb+2*sgb;
i3=[i1;i2];
i3=i3(:).';
i4=ones(1,2*sgb+2*sgb);
i5=[i4;i3];

l1=kron((9+4*sgb+1:1:9+4*sgb+2*sgb), ones(1,2));
l2=8*ones(1,2*sgb+2*sgb);
l3=[l1;l2];

m1=9+4*sgb+1:1:9+4*sgb+2*sgb;
m2=8*ones(1,2*sgb);
m3=[m1;m2];
m3=m3(:).';
m4=8*ones(1,2*sgb+2*sgb);
m5=[m4;m3];

n1=ones(1,2*sgb);
n2=10:1:9+2*sgb;
n3=[n1;n2];
n3=n3(:).';
n4=ones(1,2*sgb+2*sgb);
n5=[n4;n3];



fillAeq=[1 -1 -1 ... %observational equivalence
         1 -1    ...
         1       ...
         ones(1,numzeros) ... %zeros
         1 ...%ones 
         repmat([1 1],1, (3+sgb*3)*3) ...%symmetry
         ones(1,3)...
         repmat([1 -1], 1, (2*7)+(2*2*sgb)*3) ]; %identical      
end